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Abstract 

We show how information on the uniformity properties of a point set em- 
ployed in numerical multidimensional integration can be used to improve the 
error estimate over the usual Monte Carlo one. We introduce a new mea- 
sure of (non-) uniformity for point sets, and derive explicit expressions for the 
various entities that enter in such an improved error estimate. The use of 
Feynman diagrams provides a transparent and straightforward way to com- 
pute this improved error estimate. 
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1 Introduction 



In numerical integration, an object of paramount importance is the integra- 
tion error, that is, the difference between the actual (unknown) value of the 
integral and the obtained numerical estimate. Since exact knowledge of the 
error would imply exact knowledge of the integral (and hence would obviate 
the need for a numerical estimate), in practice one has to settle for an esti- 
mate of this error, which must be shored up by plausibility arguments. In 
deterministic integration, where the set X^r of N integration points is deter- 
mined beforehand, a deterministic error bound (that is, a true upper limit on 
the error) can often be obtained, provided the integrand falls in some known 
class of functions with, for instance, given smoothness properties. The inte- 
gration rule {i.e. the point set X^) can then be optimized to guarantee a 
rapid decrease of the error with increasing [|I|. 

Unfortunately, in many practical problems the integrand is not particu- 
larly smooth. This happens, for instance, in particle phenomenology, where 
integrands may contain discontinuities corresponding to experimental cuts 
in the admissible phase space. In such cases, one usually relies on Monte 
Carlo integration, where the point set Xjy is assumed to be chosen from the 
ensemble of random point setsQ. In the derivation of the error estimate, we 
use the fact that the point set Xn is a 'typical' member of the ensemble 
of all such point sets, and average over this ensemble. The error estimate 
is, then, probabilistic rather than deterministic; but with some care reliable 
estimates can be obtained, as implied by the various laws of large numbers. 
The Monte Carlo method can be applied to the very large class of square 
integrable integrands; on the other hand, there is the drawback that the er- 
ror decreases asymptotically only as l/y/N. The existing variance-reducing 
techniques can at most decrease the coefficient of the error estimate, and not 
its behaviour with A^, which is a direct consequence of the random structure 
of the point set. 

Recently, the approach known as Quasi-Monte Carlo has received con- 
siderable attention. Here, one attempts to construct points sets Xn that 
are 'more uniform' (in terms of a given measure of uniformity) than truly 
random ones. The existence of error bounds a la Koksma-Hlawka un- 
derlies the belief that such quasi-random point sets can, indeed, lead to an 
error that decreases faster than 1/\/N: the Roth bound suggests that a 
behaviour as (logA^)'^/A^ may be possible, where c is some constant depend- 

"'^In practice, the points are pseudo-random rather than truly random, being generated 
by a deterministic algorithm: however, we shall assume that the random number generator 
employed is sufficiently well behaved to satisfactorily mimic a truly random sequence. 
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ing on the dimensionality of the integral Many quasi-random sequences 
with asymptotically small non-uniformity have been proposed, such as van 
der Corput/Halton [Q, Faure 0, Sobol'l^, and Niederreiter ^ sequences. 
These generally perform to satisfaction in many applications: but, and this is 
the central problem addressed in this and subsequent papers, no reasonable 
way is yet known to estimate the error while making use of the improved 
uniformity of the point set. In fact, it is customary to just compute the error 
estimate as if the point set X^r were truly random. This procedure tends to 
overestimate the error, as shown by various case studies [§ , so it is certainly 
safe to do so: but, in our opinion, Quasi-Monte Carlo methods will come into 
their own only when improved error estimates are available. It is our pur- 
pose to provide an attempt at such an improved integration error estimate. 
To do this, we shall have to face the fact that quasi-random point sets are 
not 'typical' members of an obvious ensemble: indeed, they are very special, 
with a lot of ingenuity going in their construction. The central idea of our 
approach is the definition of a suitable ensemble of quasi-random point sets. 



The lay-out of this paper is as follows. In section 2, we shall discuss the 
error estimate for Monte Carlo, in such a way that it becomes clear precisely 
which information on the quasi-random point set is needed to improve the 
error estimate. To this end, we shall have to employ some definition of 
(no n-) uniformity, that is, a discrepancy. In section 3, we shall introduce 
a particular definition of a discrepancy, which we feel is better suited to 
this kind of problem than the best-known discrepancy, the so-called star 
discrepancy D^, discussed extensively in the literature [|12]. In section 4, we 
shall discuss in detail how the various ingredients in a new error estimate can 
be computed, for a point set with a given, known, value of the discrepancy. 
We shall do this by essentially combinatorial arguments, using a technique 
based on Feynman diagrams. It will appear that, in the limit of large A^, 
our approach is closely related to the computation of path integrals in field 
theory. 

Finally, we have to apply in practice what we have learned in theory. 
We shall make such attempts, in the academic but instructive case of one- 
dimensional integral, and in more dimensions: these points will be treated 
elsewhere B. 
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2 Estimates, ensembles and discrepancy 



2.1 The improved error estimate 

To set the stage, we define the D-dimensional integration region K to be 
[0, 1)^, the archetypical unit hypercube. The integrand is a function f{x) 
defined on K, which may contain discontinuities but must be at least square 
integrable. Its integral is denoted by Ji, where 

jp = j dx {f{x)y , p = i,2,... . (1) 

K 

The point set consists of D-dimensional points x^, k — 1,2,. . . ,N. 
Where necessary, we shall denote individual components of points by Greek 
indices, x^, = 1,2, ... ,D. The numerical estimate of the integral is then 
given by 

^ = ^ E fM > (2) 

k=l 

and the error made is then simply 7] = S — J\. The salient fact about Monte 
Carlo is that the point set Xjv is a random object, and so, therefore, is the 
error r}. The standard Monte Carlo estimate is derived by assuming that Xn 
is a 'typical' member of the ensemble of point sets, governed by a probability 
density Pn{xi, X2, . . . ,xn)- We shall also define marginal densities: 

Pk{xi,X2,...,Xk) = j dxk+idxk+2---dxN Pn{xi,X2,...,xn) , (3) 

K 

for A; = 0, 1, 2, ... , N, so that Pq = 1. For truly random points, we have the 
ideal iid uniform distribution: 

Pn{xi,X2,...,xn) = 1 ■ (4) 
We want, however, to be more general, and we shall write 

Pkixi,X2,...,Xk)^l--^Fk{xi,X2,...,Xk) . (5) 

Obviously, none of the Fk can exceed A^, and for truly random points they 
are identically zero. Moreover, since the order in which the points enter is 
immaterial in this integration problem, we shall assume that all Pk and 
are invariant under any permutation of the arguments. 
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We are now ready to estimate the error, by computing the various mo- 
ments of the probabihty distribution of r). Denoting by brackets the average 
over Pjv, we have for the first moment 

1 ^ 1 r 

iv) = -^i + ^E(/(^fc)) = ^ dxf{x)F,{x) . (6) 
k=i ^ 

We see that, if the integration is to be unbiased for all integrands, we needs 
must have 

Fi(x) = , (7) 

which means that, individually, each integration point Xk in must be 
uniformly distributed, and the difference between quasi-random and truly 
random point sets may show up only in correlations between points. Assum- 
ing this to be indeed the case, we then turn to the second moment: 

2 T 1 ^ 

{ri') = Jl-^T.U{x,)) + —Y.{fMf{xi)) 

k=i k.i=i 

= ^ - Jl - (l - ^) / ^^1^^2 /(xi)/(x2)F2(xi,X2)j . (8) 

The following facts become apparent from this result. In the first place, for 
truly random points, F2 vanishes and we recover the standard Monte Carlo 
estimate. Secondly, only a small, 0{1/N), deviation in Pjv from the truly 
random uniform iid case can already significantly diminish the error, due 
to the delicate cancellations involved in the computation of {rf). Finally, 
the mechanics behind this improvement become obvious: in very uniform 
point sets (such as quasi-random point sets), that have a low discrepancy, 
the points are spread 'more evenly' than in typical random sets: they 'repel' 
each other, leading to a positive value of ^2(3^1, 0:2) whenever xi and X2 are 
'close' in K: this suppresses the contribution to {rf') from regions where 
f{xY tends to be large. Conversely, for very non-uniform sets, where the 
points are more 'clumped' together, and whose discrepancy is large, F2 will 
be negative for neighbouring points, with a corresponding punishment in the 
error estimate. A very simple and instructive illustration of this idea is given 
in Appendix B. 

A final remark is in order here: for truly random point sets, the usefulness 
of the error estimate relies on the fact that the distribution of 77 tends to a 
Gaussian, as implied by the Central Limit Theorem. In principle, we ought 
to reconstruct a proof of this theorem for the more general Pjv discussed here. 
Although we have not yet done so, we have good reason to believe that a 
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Central Limit Theorem holds in our case as well, with an implied Gaussian 
width following from Eq.(|D. 



2.2 Generating functions 

Our task is now clear: we have to find, for quasi-random point sets, a workable 
definition of Pn, and a corresponding formula for F2. To do so, we first 
assume that there exists some measure of (no n-) uniformity of point sets X^: 
that is, there must be given some discrepancy a.s a function of the Xk, 
and we shall also assume that we can, for a given X^v, compute its value 
which we shall denote by s: 



A 'small' value of s shall indicate a point set that is relatively uniform com- 
pared to a truly random one. We shall defer an explicit definition of Dn 
to the next section, and for the moment just assume that one is given. We 
then propose to use for Pjy the probability density obtained by assuming that 
all points are uniformly distributed, with the additional constraint that the 
discrepancy takes on the value s: 



Hk{s;xi, . . . ,Xk) = dxk+i---dxN 5{Dn{xi,...,xn) - s) , (10) 



where we have introduced the Dirac delta distribution to handle the discrep- 
ancy constraint. The function Hq{s) is, then, the probability distribution of 
s over the ensemble of all truly random point sets, which is an interesting 
quantity in its own right. The idea behind this definition is the following. 
We are not allowed, in principle, to consider a quasi-random point set as a 
'typical' one in the whole random ensemble, since s is (hopefully) small com- 
pared to the expected value for random sets: but, in the subset with that 
particular value of s, it may very well be typical (indeed, this is the same 
fingers- crossed attitude that allows us to use, for pseudo-random point sets, 
the standard Monte Carlo error estimate in the first place). Moreover, if no 
information whatsoever is available on the value of s, we have to integrate 
over all s with probability density Hq{s), upon which the delta function con- 
straint drops out and we are back in the truly random case. 



Dn{xi,X2, . . .,Xn) = s 



(9) 



PAr(s;Xi, ...,Xn) 



Hn{s] Xi, . . .,Xn) 

Hois) 



K 
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We now proceed to calculate F2. Let us define a set of moment- generating 
functions as follows: 



Gk{z;xi,...,Xk) = j dxk+i---dxN ew{zDN{xi,...,XN)) , (H) 

K 

that is, the generating function where the first k of the integration points 
are kept fixed, and the remaining ones are integrated over: Gq{z), then, is 
the moment-generating function for Ho{s). We employ the definition of the 
Dirac delta distribution as a Fourier/ Laplace transform to write 

Hk{s;xi,...,Xk) = — j dz e~^^Gk{z;xi,...,Xk) , (12) 

—ioo 

where the integration contour must run to the left of any singularities. So, 
knowledge of the allows us to compute everything we need. We may 
uniquely split each G^ into a constant part, and a part that depends on the 
Xi,...,k and averages out to zero: 

Gk{z]Xi,...,Xk) = Go{z) + ^pk{z;xi,. . . ,Xk) , 

j dxi- ■ -dxk Pk{z;xi, . . . ,Xk) = 0. (13) 

K 

It follows that we can express as 

Fk{s;xi,...,Xk) = -Rkis;xi,...,Xk)/Ho{s) , 

^ ioo 

Rk{s;xi,...,Xk) = ^ / dz e~^^pk{z;xi,...,Xk) . (14) 



-too 



In this way, we can compute the necessary two-point correlation F2, provided 
we can compute, at least in the form of an asymptotic expansion in 1/A^, the 
moment-generating functions G^'- the dominant contribution to will then 
come from the leading terms in Gq and pk- This depends, of course, on a 
manageable choice of the discrepancy Dn, and this will be the subject of the 
next section. We wish to point out that, actually, an exact expression for 
the F|^. would allow us to determine the lowest possible value for s, simply 
by putting F^^s; Xi, . . . , Xk) = N; from the fact that it is so hard to improve 
on the Roth bound, we may infer that an asymptotic expansion in 1/A^ is 
probably the best possible result we can hope for at this moment. Finally, 
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we remark that a check on our results for R2{s;xi,X2) is provided by the 
identities 

oo 

J dx2 R2is;xi,X2) = , J ds R2{s; xi, X2) = , (15) 

K 

where the first identity follows from Eq. (|^) , and the second one by construc- 
tion. 



3 The Fourier discrepancy 

We now turn to a useful definition of a discrepancy measure for X^r. This 
we do by first constructing a model of the class of integrands likely to be 
encountered, and then considering the squared integration error expected for 
this class of integrands. This procedure is analogous to that employed by 
Wozniakowski [p!0[] . 



3.1 The one-dimensional case 

For simplicity, we shall first discuss the case D = 1, and only further on 
generalize to more dimensions. We start by assuming that our integrands 
admit of a decomposition into a set of orthonormal functions. These satisfy 

dx U^{x)Un{x) = 6mn , ^n{x)Un{y) = S{x - y) . (16) 

A general integrand f{x) can then be written as 

n>0 

SO that the coefficients f„ determine the function. Since these coefficients 
form an enumerable set, it is fairly easy to set up a combined probability 
distribution for them, which gives us then a probability measure Vf on the 
class of integrands. We take the following choice: 



so that each coefficient f„ is distributed normally with mean value and 
standard deviation cr„. We shall call cr„ the strength of the mode Un{x) in 
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our function ensemble. Denoting by () / an average over the ensemble under 
the measure Vf, we then have 

{Vn)f = , (VnVm) f = al^Smn ■ (19) 

We can also compute the moments of the integration error 77 over the inte- 
grand ensemble. The integral Ji is of course just given by Vq, and moreover 
we have 

1 ^ 

k,l=l 

where we have introduced functions Pm as 

Pm{x,y) = Yl '^'rrUn{x)uM ■ (21) 

These functions, which will play an important role in our discussion, are 
symmetric in their arguments, and have the properties that 

j dy Prn{x,y) ^0 , j dy Prn{xi,y)Pn{x2,y) ^ Pm+n{xi,X2) , (22) 
K K 

for m, n > 0. In fact it is also simple to prove that the error rj has a Gaussian 
distribution, not as a consequence of any law of large numbers, but because 
of the Gaussian character of our measure Vf on the space of integrands. 

Up to this point, we have not had to specify the particular orthonormal 
functions. In practice, we shall use the Fourier basis: 

uo{x) = 1 , U2n-i{x) = V2 sin(27ma;) , U2n{x) = V2 cos(27rna;) , (23) 

where n runs over the positive integers. 

In addition, we shall make one very important assumption, namely that 
the sines and cosines have equal strength: 

(^2n-l = Cr2n ■ (24) 

We shall call this property translational invariance. Its physically reasonable 
motivation is the fact that under this assumption the mode with frequency 
n (made up from U2n-i and U2n) has a uniformly distributed phase. Trans- 
lational invariance will enable us to considerably simplify our results. Most 
importantly, it leads us to write 

(3m(xi,X2) = Y.'^((^2n?"'cos{27rn{xi-X2)) , (25) 

n>0 
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so that the functions (3 only depend on the difference between xi and X2, and 
we may also write f3m{xi — ^2). Another attractive consequence is the fact 
that (as we shall see) point sets that only differ by a translation (assuming a 
periodic extension of K and Xjv) have the same discrepancy, in contrast to, 
e.g., the star discrepancy D^^. 

The above consideration leads us to propose, as an appropriate definition 
of discrepancy, the following: 

1 ^ 

Dn = Dn{xi,. . . ,xn) = — Pi{xk,xi) . (26) 

k,l=l 

Note that we have taken out one factor of this will make the discrepancy 
independent of for truly random points. It is easily seen that, for truly 
random points, 

(Dn) = E • (27) 

n>0 

There are a few additional observations in order here. In the first place, 
note that the roles of integrand and point set have, in some sense, been 
reversed here: in the derivation of the standard Monte Carlo error, we keep 
the integrand fixed, and average over the ensemble of point sets, upon which 
(?7^) no longer depends on the particular point set but only upon a property 
of the integrand (namely its variance); while here we have kept the point set 
fixed, and averaged over an ensemble of integrands, so that (77^) / no longer 
depends on the integrand, but only upon a property of the point set Xjy, 
namely its discrepancy. A number of other results along these lines have 
been presented in [^. In the second place, the zero mode with n = does 
not enter in the discrepancy: this is reasonable since vq is in fact the desired 
integral, and the error just consists in our partial inability to get rid of the 
higher modes, with n > 0. Finally, the strengths ct„ cannot be chosen too 
arbitrarily. If our average integrand is to be quadratically integrable, we need 
to have 

= E^n<oo , (28) 

/ 

while additional smoothness requirements will ask for an even faster conver- 
gence of this sum. An admissible, and reasonable, choice for the (T„ will be, 
for instance 

(70 = 1 , cr2n = l/n , n>0 , (29) 
and this we shall consider in the practical applications P]. 
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3.2 The more- dimensional case 



The extension of the above discussion to more dimensions is quite straight- 
forward. We shall use again x to denote £)-dimensional points. The complete 
set of orthonormal functions is now enumerated not by a scalar n, but by 
a vector quantity n with indices > 0, with // = 1,2, . . . , D. The set of 
orthonormal functions is now given by 

D 

Unix) = Hun.ix'') ; (30) 

the measure Vf is straightforward, with the strengths of the various modes 
denoted by a^, and the discrepancy is again given by 

1 ^ 

Dn^ Dn{xi,...,xn) = Yl Pii^k,xi) , (31) 

k,l=l 

where 

Pm{x,y)^J2^'frMx)'^n{y) . (32) 

n>0 

Here, the sum runs over all n except n = (0, 0, ...,0). Again, for truly 
random points we have 

{Dn) = E4 ■ (33) 

n>0 

Translational invariance in D dimensions requires 

<^(2ni,2n2,...,2n^) — Cr(2ni-l,2n2,...,2n^) = " ' ' 
■ ■ ■ — '^(2n^,2n^,...,2nO-l) = ' ' ■ 

••• = 0-(2ni-l,2n2-l,...,2nB-l) , (34) 

so that the strengths must be equal in groups of 2^^'^ , where k is the number 
of vanishing components of n. Again, for our integrands to be quadratically 
integrable on the average, we must have 

E4<~ . (35) 

n>0 

Now, we may choose to let the strengths be dominated by the highest partial 
frequency, for instance 

— c(maxn'') . (36) 

Because of the increasing multiplicity with max^n'*, cr^j must then decrease 
more rapidly than in the one-dimensional case. Other reasonable choices. 
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such as an = crij^^n'^), lead to the same conclusion. As usual, we encounter 
here the phenomenon that a smooth function in one dimension is not pre- 
cisely the same concept as the projection onto one dimension of a function 
that is smooth in many dimensions: the 'curse of dimensionality' crops up 
again, in disguise. 



We have now finished the first part of our program, that is, the establish- 
ing of a reasonable definition of a discrepancy. We want to remark that it 
differs from the better known star- discrepancy D^, discussed so extensively 
r2|. That one can, in fact, also be derived from a class of integrands. 



m 



as first made explicit in |]I0|. The appropriate integrand class is defined by 
taking for Vf the Wiener sheet measure; although mathematically attrac- 
tive, this is by no means the preferred choice for many practical applications, 
since it singles out functions that are everywhere continuous but nowhere 
different iable. Moreover, the translational invariance will lead, as we shall 
see, to Fi{s; xi) = 0, and the lack of this invariance for the star-discrepancy 
probably implies that we cannot simply prove that it leads to unbiased inte- 
gral estimates. We should like to point out that, in fact, we have been able 
to derive the form of Hq{s) for the star-discrepancy as well |p!3 |. 



4 Discrepancies by Feynman diagrams 

We must now proceed with our program, and find expressions for 
Gkiz;xi,...,Xk)=22 — 7{^N)k , 

1 f ^ 

J dXk+l---dXM ^ Xra) ■ ■ ■ /5i(Xr2„_i, Xra^) . (37) 



^•l,...,2m = l 

We shall call the points Xi, . . . ,Xk, that are kept fixed, external points, and 
the remaining N — k ones, that are integrated over, internal points. Note 
that, in the multiple sum, each index runs over both external and internal 
points. The various internal points are essentially indistinguishable, and they 
give rise to combinatorial factors; the study of these factors will enable us to 
write (-Dtv)^ an asymptotic series in 

The part of the multiple sum where all indices ri , . . . , r2m are equal and 
internal gives rise to a single combinatorial factor N—k; if they are distributed 
over two different internal values, we shall have a combinatorial factor {N — 
k) {N — k — 1), and so on. It follows that the largest combinatorial factor that 
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can, in principle, occur, is the one for the case where all indices are internal 
and distinct. However, that contribution consists of separate factors 



dxry-^ dxy. I3\ {x^^ , x^. ) 



K 



which vanish. Therefore, the actually largest combinatorial factor is that for 
the situation where the indices are all internal, and fall in distinct pairs: the 
corresponding combinatorial factor is 



{N-k)'^ = {N -k)\/{N -k-m)\ 



m— 1 



2 

m- 



h km 

2 



The 'falling power' notation is taken from Graham et al. |jTj]; its asymptotic 
expansion for large is derived in Appendix A. 

4.1 Feynman diagrams 

We are now ready to compute the first few terms in the expansion of 
{D^)k- We start by introducing a diagrammatical technique. Each internal 
point we shall denote by a dot, and each external point by a cross. Every 
function /3i shall be denoted by a link, a solid line between dots or crosses. 
Every contribution to (-D]v)fc ^ill therefore contain precisely m links. Identi- 
cal values of the summation indices mean that the corresponding points will 
be contracted, leading to vertices with two, three, four, . . . legs. Neglecting 
for the moment the combinatorial we have, for instance, 

^ = Pl{Xi,X2) , 

Xl X2 

dy /3i{x,y) = , 
dyidy2 (3i{yi,y2) = , 

dy f3i{y,y) , 



X 



K 



K 



o 



K 



— = J dyidy2 f3i{yi,yi)l3i{yi,y2)l3i{y2,y2) , (39) 

K 

where the two zero results follow from Eg. (^2]) . In each diagram or product 
of diagrams, the combinatorial factor can be read off immediately from the 
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number of internal points in evidence: if there are p internal points, the factor 
is (A^ — k)^. External points do not contribute a combinatorial factor, but it 
must be kept in mind that we shall have to sum over all external points. For 
instance, the diagram x — x must be interpreted as 

N 

>^ = ^ X X , (40) 

fc,l=i -^k 
k^l 

so that this diagram actually contains k{k — 1) terms. 

In writing out the expansion, we can simply determine which kinds 
of diagrams will contribute, as follows. The leading power of a two-point 
internal vertex is 1 , that of a three-point internal vertex is 1 / \/iV, for a four- 
point internal vertex we have and so on. Formally, this is similar to 
a field theory with a universal coupling constant g proportional to 1/ -\/iV. 
Moreover, each external point will effectively carry a factor 1 / ^/N. If we de- 
cide to keep only the terms up to and including we shall have to allow 
for at most two external points, or a single external point and one internal 
three-vertex, or two internal three-vertices, or a single internal four-vertex. 
In a field theory, this corresponds to the first-order correction, proportional 
to g"^. Moreover, our property of translational invariance implies precisely 
the analogue for a field theory, namely momentum conservation^. Also, the 
external points, that act somewhat like sources, need rescaling by factors 
V^, which could be envisioned as the truncation process by which a Green's 
function is turned into an S'-matrix element. This is precisely what the ex- 
traction of the factor in Eq.(^ does for us, and so we recognize what 
the function F2{s;xi,X2) actually is: it is the full propagator. 



4.2 An illustration: the first two moments 

We shall now show how to apply the diagrammatic techniques in the calcula- 
tion of Gk{z; xi, . . . , Xk). The first nontrivial term comes from {DN)k, which 
we can write diagrammatically as 



1 

N 



: + (^+2{N -k) 



^The analogy with a real field theory cannot be carried too far, however, since such a 
theory has an infinite number of degrees of freedom. In the present case, all vertices and 
combinatorial factors also carry subleading contributions, which must be properly taken 
into account. 
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+ {N -k)(^+{N -k)^ 



(41) 



We may considerably simplify this. In the first place, we have x — • = and 
• — • = 0, see Eg. (pUj) . In addition, the property of translational invariance 
has the important consequence that a diagram with only a single external 
point evaluates to a constant. For instance. 



= A(x,x) =/3i(0) = J dy(3,{y,y) = Q . (42) 

K 



This also implies that, whenever two parts of a diagram are connected by a 
single vertex, we may split it up (of course, without changing the combina- 
torial factor!), so that, for instance. 



2 



ChO - (o) 



; (43) 



the last line is an example of the more general phenomenon that all tadpole 
diagrams with only internal points evaluate to zero. Again, this is due to our 
assumption of translational invariance: incidentally, it immediately proves 
that 

Gi{z; xi) = Go{z) to all orders in 1/A^, (44) 
which in its turn implies that 

Fi{s; xi) = to all orders in 1/N, (45) 

as required in Eq.(0). We can now write {Dj^)k as 

where we must keep in mind that the second term implies a summation over 
the k{k — 1) pairs of different points in xi, X2, . . . , x^. It follows trivially that, 
as required, 

dxk {DN)k = {DN)k-i =^ 

=^ j dxk Fk{s;xi, . . . ,Xk) = Fk-i{s;xi, . . . ,Xk-i) ■ (47) 

K 
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K 



We now proceed to the next order, and compute {D]^)k. The only contribu- 
tions that do not immediately vanish under Eg .(^21) are 



1 

iV2 



+ 4 



+ 2(iV-A;)0 + + 4 

+ no + (N-k) 



(48) 



Using translational invariance, we may rewrite this as 



+ 4(A^ - A;) >^-.-^ + 4>^^ 



+ 2 



(49) 



It can again easily be checked that 

j dxk {Dl)k = (Dl) 



fc-i 



(50) 



If we restrict ourselves to terms of order O (1) and O {N ^), we have 



(51) 



4.3 The general result to leading order 

We shall now discuss the computation of all diagrams of the leading necessary 
order. To start, let us concentrate on the leading terms, of order O (1). These 
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are contributed by graphs containing only internal two-vertices, that is, they 
must be of the form of a product of closed loops: 

W„(ra) . £ AKp,.p..p3....)(0)'"(Or(Or-"<52) 

Pl,2,3,...>0 

where we have left out the factor N—/N"^, and imply the constraint 

m = pi + 2p2 + 3p3 -\ . (53) 

The factor A{m;pi,p2, ■ ■ ■) is governed by a recursion relation. We may go 
from m — 1 to m by either adding a single one-link loop, or by putting an 
extra link in any /c-link loop, thereby turning it into a {k + l)-link loop. The 
recursion relation therefore reads 

A{m;pi,p2,P3, ■■■) = - - 1,P2,P3, • • •) + 

^2k{pk + l)A{m-l;pi,p2,...,Pk + l,Pk+i-'^,---) ■ (54) 

A;>0 

We may use the Ansatz 

A{m;pi,p2,P3,. . .) = . ' . Ca^'afaf--- , (55) 

P1IP2W ■ ■ ■ 

with C and the a's to be determined. Putting this Ansatz in Eq.(|5^) leads 
us to 

m= — + Y^ 2k pk+i . (56) 

^1 fc>0 ^k+i 

Since the coefficients of the p's are known from Eg . (|53|) , this gives us the 
symmetry factor associated with each closed fc-link loop: 

2^ 

ak = — ■ (57) 
2k ^ ' 

By inspection of {Di^i)^, we also find C = 1. The result for the terms that 
give the leading contribution is therefore 

(58) 

We now come to the subleading terms. Let us introduce the notation 
12 3 



X X — X X , X X — X • X , X X 



(59) 
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so as to indicate the number of links between two vertices. The sublead- 
ing contributions are characterised by the additional presence of one of the 
following diagrams: 



oo 




It is actually only the first of these graphs that is needed for the purpose of 
computing pk to first order: the other two will only contribute to subleading 
terms in Hq{s) which are not relevant to the order we are working in here. 
The relevant contribution is, therefore 

Pk{z]Xi,. . . ,Xk) = Wi{m) + O (^^^ , 
W,{m) ^ Y: B{m;q,p„p„...)>J^{Qy'{0)''... ,(60) 



P1.2....>0 

g>0 



with 



m = q + Pi + 2p2 + 3p3 -\ . (61) 

The coefficient B satisfies a recursion relation similar to that of A: 

B{m;q;pi,p2,...) = 2qB{m- l;g- l,pi,p2,---) 
+ 5(m- l;g,pi - l,p2,---) 

+ ^ 2/i;(pfc + 1)5("^ - 1; g,pi,p2, • • • ,Pfc + l,Pfc+i - 1, • • •) • (62) 

fc>0 

By the same trick as above, we may solve this relation. The symmetry factors 
for the loop diagrams are of course the same, and the diagram > < x has a 
symmetry factor 2'^~^. We therefore find 

"■--.,F.^(?-^)(tO)"(?0)"-.» 

In the above derivations, it may be noted that the various powers of 2 arise 
from the fact that each link has two endpoints, while the remaining symmetry 
factors are just those of diagrams with topologically indistinguishable points 
(l/2n for an n-link loop, 1/2 for the propagator diagram). 



4.4 The generating functions 

We can now write down immediately the leading form for Gq and pk-, by 
taking the appropriate sums over m. To this leading order, we do not have 
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to worry about subleading terms in (A^ — k)— and {N — /c)^^^. We have 

m>U m>U 



= exp 

and 



~m 1 / \ 

Pfc(z;a;i,...,a;,) ~ ^ ^^^(rn) = ^ (2;^)^>^ Go(z) . (65) 

m>l ^ g>0 V / 

To evaluate the various diagrams we employ the functions (3. From the 
following propagator diagrams, with endpoints x and y: 

1 



2 

X X 



3 

X X 



dz (3i{x,z)(3i{z,y) = (32{x,y) 
dz /32ix, z)/3i{z, y) = p^ix, y) 



K 



K 



= Pq{x,y) = E(^")^'' ' (66) 

n>0 

we find the following representations: 

^iz;x,y) ^ Y.i'^zy^ = E . 2 ^^nix)uniy) , (67) 

9>1 n>0 ^ '^^^n 

logGo(^) = -iElog(l-2;^a^.) , (68) 

^ n>0 

and ^ 

pk{z]Xi,. . .,Xk) = -'^(t){z\Xi,Xj)GQ{z) , (69) 

with the indices i and j running from 1 to k. For the case /c = 2 in which we 
are interested, this specializes to 

p2{z;xi,X2) = (t){z;xi- X2)Gq{z) . (70) 

In the language of field theory, the function P2{z; ,xi,X2)/N may be recog- 
nized as the Dyson-summed, or dressed, two-point Green's function, including 
vacuum diagrams Go{z). Here Gq{z) plays the role of the path integral in 
the free-field approximation. 
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5 Conclusions 



In this first paper we have addressed the question of how, given informa- 
tion on the uniformity property of a particular point set in terms of a 
discrepancy (defined in Eq.(pID), we may hope to improve on the error es- 
timate when we apply the point set Xn in numerical integration. To this 
end, we need the function F2{s]Xi,X2) which is defined as the ratio of two 
functions, R2{s;xi,X2) and Ho{s), which are themselves again given as the 
Fourier transforms of pk{z;xi,X2) and Go{z), respectively. We have devel- 
oped a diagrammatic approach that allows us to systematically compute a 
series expansion for these two objects in powers of We have derived 

explicit expressions for the leading term in these expansions: they are given 
by Eq.(|5^ and Eq. (pBD , respectively. 

Of course, the above discussion has been purely formal. We have been 
able to point out interesting parallels between our results and properties of 
a classical field theory; but all this will be moot unless we can, in fact, turn 
what we have learned into an expression for F2{s] xi, X2) that is explicit and 
simple to evaluate. This will be the subject of further publications in this 
series. 

Apart from the explicit expressions for {D^)^ and (D^)^, our results 
have been strictly limited to the leading terms, that are independent of A^. 
Since any numerical integration usually employs a large value for anyway, 
we feel that we are justified in this. However, there are issues for which a 
further expansion ought to be useful. For instance, it would be interesting 
to see how far we could approximate a lower bound on D^r a la Roth, by 
either computing further terms in F2 and putting its value to A^, or by trying 
to establish bounds on that value for s where Ho{s) vanishes, for finite A^. 
These issues are beyond our scope at this moment, but we feel that we have 
at least laid some of the groundwork in this paper. 



Appendix A: Expansion of falling powers 

In this appendix we give a simple derivation of the asymptotic expansion of 
the quantity (A^ — k)—. It is based upon the binomial-theorem fact 

y —a^x"" = (1 + xY . (71) 



We therefore have 
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Now, one the one hand, we have 



1 + ^ 



N 



(73) 



and, on the other hand. 



X 



kx k{k + 



(74) 



Multiplying the expansions ( [74D and (|73D , and carefully keeping track of the 
coefficient of x"^, we then find 



(AT - k)-^ 

]\fm 



1 



2 



4 , 1^3 1 



1 

and, of course, also the subleading expressions 

1 



m- H — m- H — /cm- H — A;(A; + l)m- 
3 2 2 



(75) 



(AT - k)^/N"' 
(AT - kj^/N"^ 



1 1 



m-+ {k- l)(m - 1) 



(76) 



Obviously, this expansion can be continued ad nauseam; moreover, it also 
provides a way of computing complicated constraint sums of products like 



kik2 ■ ■ ■ k 

ki,k2,...,kp=l 



with the constraints that all indices must be different. 



Appendix B: an illustrative model 

In this appendix we show how error improvement arises in a very simple 
one-dimensional model. We assume that the functions to be integrated have 
only modes up to n, and that those modes have equal strength. That is, we 
take 

a^ = l/2 , A; = l,2,...,2n , 

al = Q , k>2n . (77) 



20 



For this extremely simple case, we have 

1 2^ " 

Go{z) = , (f){z;x)^-—^cos{2nkx) . (78) 

We may easily compute the Laplace/Fourier transforms, and find 

1 / s\ " 

Ho{s) ^ -——s^-'e-' ' ^2(s;x) = 2 1-- 5]cos(27rA;x) . (79) 
(77; i). \ ny 



The average discrepancy is (s) = n for truly random points. A typical 
integrand is, in this model, of the form 

n 

fix) = J2^k cos(27rfc(x + ak)) , (80) 

k=l 

with arbitrary Vk and ak- We immediately find 

1 " 

Ji = , J2^^Y.^l ' (81) 

and 

/ dxidx2 /(xi)/(x2)F2(s; -X2)^(l--)lf2vl . (82) 

We conclude that, with the inclusion of the two-point correlation function F2, 
the straightforward Monte Carlo error estimate is changed into the Quasi- 
Monte Carlo one, as follows: 

(?7^) Monte Carlo — 

— (?7^) Quasi-Monte Carlo = ■^(?7^) Monte Carlo + ^ (^J^p^ ■ (^3) 

We see that the error is improved if the actual discrepancy is small compared 
to its expected value. 
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